/* File: fem_adapt_algs.C
 * Authors: Terry Wilmarth, Nilesh Choudhury
 * 
 */


/** This module implements high level mesh adaptivity algorithms that make use 
    of the primitive mesh adaptivity operations provided by fem_adapt(_new).
*/

#include "ParFUM.h"
#include "ParFUM_internals.h"


#define MINAREA 1.0e-18
#define MAXAREA 1.0e12

#define GRADATION 1.2
//#define ADAPT_VERBOSE

CtvDeclare(FEM_Adapt_Algs *, _adaptAlgs);

FEM_Adapt_Algs::FEM_Adapt_Algs(FEM_Mesh *m, femMeshModify *fm) 
{ 
  theMesh = m; 
  theMod = fm; 
  //theAdaptor = theMod->fmAdapt;
  theAdaptor = theMod->fmAdaptL;
}
FEM_Adapt_Algs::FEM_Adapt_Algs(femMeshModify *fm) {
  theMesh = NULL;
  theMod = fm;
  theAdaptor = theMod->fmAdaptL;
}

FEM_Adapt_Algs::~FEM_Adapt_Algs() {
}



/** Perform refinements on a mesh.  Tries to maintain/improve element quality
 * as specified by a quality measure qm;
 * if method = 0, refine areas with size larger than factor down to factor
 * if method = 1, refine elements down to sizes specified in sizes array
 * Negative entries in size array indicate no refinement. 
 */
void FEM_Adapt_Algs::FEM_Refine(int qm, int method, double factor, 
				double *sizes)
{
  SetMeshSize(method, factor, sizes);
  GradateMesh(GRADATION);
  (void)Refine(qm, method, factor, sizes);
}

/** The actual refine in the previous operation
 */
int FEM_Adapt_Algs::Refine(int qm, int method, double factor, double *sizes)
{
  // loop through elemsToRefine
  int elId, mods=0, iter_mods=1;
  int elemWidth = theMesh->elem[0].getConn().width();
  refineElements = refineStack = NULL;
  refineTop = refineHeapSize = 0;
  int elemConn[3];
  int count = 0;
  while (iter_mods != 0 && count<20) {
    count++;
    iter_mods=0;
    numNodes = theMesh->node.size();
    numElements = theMesh->elem[0].size();
    // sort elements to be refined by quality into elemsToRefine
    if (refineStack) delete [] refineStack;
    refineStack = new elemHeap[numElements];
    if (refineElements) delete [] refineElements;
    refineElements = new elemHeap[numElements+1];
    for (int i=0; i<numElements; i++) { 
      if (theMesh->elem[0].is_valid(i)) {
	// find maxEdgeLength of i
	double tmpLen, avgEdgeLength=0.0, maxEdgeLength = 0.0;
	theMesh->e2n_getAll(i, elemConn);
	bool notclear = false;
	for(int j=0; j<elemWidth; j++) {
	  int nd = elemConn[j];
	  if(nd>=0) notclear = theMod->fmLockN[nd].haslocks();
	  //else cannot look up because connectivity might be screwed up,
	  //I am hoping that at least one of the nodes will be local
	}
	if(notclear) continue;
	for (int j=0; j<elemWidth-1; j++) {
	  for (int k=j+1; k<elemWidth; k++) {
	    tmpLen = length(elemConn[j], elemConn[k]);
	    if(tmpLen < -1.0) {notclear = true;}
	    avgEdgeLength += tmpLen;
	    if (tmpLen > maxEdgeLength) maxEdgeLength = tmpLen;
	  }
	}
	if(notclear) continue;
	avgEdgeLength /= 3.0;
	double qFactor=getAreaQuality(i);
	if (theMesh->elem[0].getMeshSizing(i) <= 0.0) 
	  CkPrintf("WARNING: mesh element %d has no sizing!\n", i);
	if ((theMesh->elem[0].getMeshSizing(i) > 0.0) &&
	    (avgEdgeLength > (theMesh->elem[0].getMeshSizing(i)*REFINE_TOL))){
	    //|| (qFactor < QUALITY_MIN)) {
	  Insert(i, qFactor*(1.0/maxEdgeLength), 0);
	}
      }
    }
    while (refineHeapSize>0 || refineTop > 0) { // loop through the elements
      int n1, n2;
      if (refineTop>0) {
	refineTop--;
	elId=refineStack[refineTop].elID;
      }
      else  elId=Delete_Min(0);
      if ((elId != -1) && (theMesh->elem[0].is_valid(elId))) {
	int n1, n2;
	double tmpLen, avgEdgeLength=0.0, maxEdgeLength = 0.0;
	int maxEdgeIdx=0;
	theMesh->e2n_getAll(elId, elemConn);
	bool notclear = false;
	for(int j=0; j<elemWidth; j++) {
	  int nd = elemConn[j];
	  if(nd>=0) notclear = theMod->fmLockN[nd].haslocks();
	  //else cannot look up because connectivity might be screwed up,
	  //I am hoping that at least one of the nodes will be local
	}
	if(notclear) continue;
	for (int j=0; j<elemWidth-1; j++) {
	  for (int k=j+1; k<elemWidth; k++) {
	    tmpLen = length(elemConn[j], elemConn[k]);
	    if(tmpLen < -1.0) {notclear = true;}
	    avgEdgeLength += tmpLen;
	    if (tmpLen > maxEdgeLength) { 
	      maxEdgeLength = tmpLen;
	      if(k==elemWidth-1) {
		maxEdgeIdx = j+1;
	      }
	      else maxEdgeIdx = j;
	      n1 = elemConn[j]; n2 = elemConn[k];
	    }
	  }
	}
	if(notclear) continue;
	avgEdgeLength /= 3.0;
	//double qFactor=getAreaQuality(elId);
	if ((theMesh->elem[0].getMeshSizing(elId) > 0.0) &&
	    (avgEdgeLength>(theMesh->elem[0].getMeshSizing(elId)*REFINE_TOL))){
	  //|| (qFactor < QUALITY_MIN)) {
	  //decide if we should do a flip or bisect
	  bool flag = flipOrBisect(elId,n1,n2,maxEdgeIdx,maxEdgeLength);
	  if(flag) {
#ifdef DEBUG_FLIP
	    if(theAdaptor->edge_flip(n1, n2) > 0)  iter_mods++;
#endif
	  }
	  else {
	    if(theAdaptor->edge_bisect(n1, n2) > 0)  iter_mods++;
	  }
	}
      }
      CthYield(); // give other chunks on the same PE a chance
    }
    mods += iter_mods;
#ifdef ADAPT_VERBOSE
    CkPrintf("[%d]ParFUM_Refine: %d modifications in last pass.\n",theMod->idx,iter_mods);
#endif
#ifdef DEBUG_QUALITY
    tests(false);
#endif
  }
#ifdef ADAPT_VERBOSE
  CkPrintf("[%d]ParFUM_Refine: %d total modifications.\n",theMod->idx,mods);
#endif
  delete[] refineStack;
  delete[] refineElements;
  return mods;
}

/** Perform coarsening on a mesh.  Tries to maintain/improve element quality
 * as specified by a quality measure qm;
 * if method = 0, coarsen areas with size smaller than factor up to factor
 * if method = 1, coarsen elements up to sizes specified in sizes array
 * Negative entries in size array indicate no coarsening. 
 */
void FEM_Adapt_Algs::FEM_Coarsen(int qm, int method, double factor, 
				 double *sizes)
{
  SetMeshSize(method, factor, sizes);
  GradateMesh(GRADATION);
  (void)Coarsen(qm, method, factor, sizes);
}

/** The actual coarsen in the previous operation
*/
int FEM_Adapt_Algs::Coarsen(int qm, int method, double factor, double *sizes)
{
  // loop through elemsToRefine
  int elId, mods=0, iter_mods=1, pass=0;
  int elemWidth = theMesh->elem[0].getConn().width();
  double qFactor;
  coarsenElements = NULL;
  coarsenHeapSize = 0;
  int elemConn[3];
  while (iter_mods != 0) {
    iter_mods=0;
    pass++;
    numNodes = theMesh->node.size();
    numElements = theMesh->elem[0].size();
    // sort elements to be refined by quality into elemsToRefine
    if (coarsenElements) delete [] coarsenElements;
    coarsenElements = new elemHeap[numElements+1];
    coarsenElements[0].len=-2.0;
    coarsenElements[0].elID=-1;
    for (int i=0; i<numElements; i++) { 
      if (theMesh->elem[0].is_valid(i)) {
	// find minEdgeLength of i
	theMesh->e2n_getAll(i, elemConn);
	bool notclear = false;
	for(int j=0; j<elemWidth; j++) {
	  int nd = elemConn[j];
	  if(nd>=0) notclear = theMod->fmLockN[nd].haslocks();
	  //else cannot look up because connectivity might be screwed up,
	  //I am hoping that at least one of the nodes will be local
	}
	if(notclear) continue;
	double tmpLen, avgEdgeLength=0.0, 
	  minEdgeLength = length(elemConn[0], elemConn[1]);
	for (int j=0; j<elemWidth-1; j++) {
	  for (int k=j+1; k<elemWidth; k++) {
	    tmpLen = length(elemConn[j], elemConn[k]);
	    if(tmpLen < -1.0) {notclear = true;}
	    avgEdgeLength += tmpLen;
	    if (tmpLen < minEdgeLength) minEdgeLength = tmpLen;
	  }
	}
	if(notclear) continue;
	avgEdgeLength /= 3.0;
	qFactor=getAreaQuality(i);
	if (((theMesh->elem[0].getMeshSizing(i) > 0.0) &&
	     (avgEdgeLength < (theMesh->elem[0].getMeshSizing(i)*COARSEN_TOL)))
	     || (qFactor < QUALITY_MIN)) {
	    //){
	  //CkPrintf("Marking elem %d for coarsening\n", i);
	  Insert(i, qFactor*minEdgeLength, 1);
	}
      }
    }
    while (coarsenHeapSize>0) { // loop through the elements
      elId=Delete_Min(1);
      if ((elId != -1) && (theMesh->elem[0].is_valid(elId))) {
	theMesh->e2n_getAll(elId, elemConn);
	bool notclear = false;
	for(int j=0; j<elemWidth; j++) {
	  int nd = elemConn[j];
	  if(nd>=0) notclear = theMod->fmLockN[nd].haslocks();
	  //else cannot look up because connectivity might be screwed up,
	  //I am hoping that at least one of the nodes will be local
	}
	if(notclear) continue;
	int n1=elemConn[0], n2=elemConn[1];
	double tmpLen, avgEdgeLength=0.0, 
	  minEdgeLength = length(n1, n2);
	for (int j=0; j<elemWidth-1; j++) {
	  for (int k=j+1; k<elemWidth; k++) {
	    tmpLen = length(elemConn[j], elemConn[k]);
	    if(tmpLen < -1.0) {notclear = true;}
	    avgEdgeLength += tmpLen;
	    if (tmpLen < minEdgeLength) {
	      minEdgeLength = tmpLen;
	      n1 = elemConn[j]; n2 = elemConn[k];
	    }
	  }
	}
	if(notclear) continue;
	CkAssert(n1!=-1 && n2!=-1);
	avgEdgeLength /= 3.0;
	qFactor=getAreaQuality(elId);
	// coarsen element's short edge
	if (((theMesh->elem[0].getMeshSizing(elId) > 0.0) &&
	     (avgEdgeLength < (theMesh->elem[0].getMeshSizing(elId)*COARSEN_TOL)))
	      || (qFactor < QUALITY_MIN)) {
	  //){

	  int eNbr = theMesh->e2e_getNbr(elId, theMesh->e2n_getIndex(elId,n1)+
					 theMesh->e2n_getIndex(elId,n2)-1, 0);
	  // determine if eNbr should also be coarsened
	  if ((eNbr >= 0) && (theMesh->elem[0].is_valid(eNbr))) {
	    theMesh->e2n_getAll(eNbr, elemConn);
	    avgEdgeLength=0.0;
	    for (int j=0; j<elemWidth-1; j++) {
	      for (int k=j+1; k<elemWidth; k++) {
		avgEdgeLength += length(elemConn[j], elemConn[k]);
	      }
	    }
	    avgEdgeLength /= 3.0;
	    qFactor=getAreaQuality(eNbr);
	    if (((theMesh->elem[0].getMeshSizing(eNbr) > 0.0) &&
		 (avgEdgeLength < (theMesh->elem[0].getMeshSizing(eNbr))))
		|| (qFactor < QUALITY_MIN)) {
	      //CkPrintf("Coarsening elem %d which has desired edge length %6.6e and average edge length %6.6e.\n", elId, theMesh->elem[0].getMeshSizing(elId), avgEdgeLength);
	      if (theAdaptor->edge_contraction(n1, n2) > 0)  iter_mods++;
	    }
	  }
	  else {
	    //CkPrintf("Coarsening elem %d which has desired edge length %6.6e and average edge length %6.6e.\n", elId, theMesh->elem[0].getMeshSizing(elId), avgEdgeLength);
	    if (theAdaptor->edge_contraction(n1, n2) > 0)  iter_mods++;
	  }
	}
      }
      CthYield(); // give other chunks on the same PE a chance
    }
    mods += iter_mods;
#ifdef ADAPT_VERBOSE
    CkPrintf("[%d]ParFUM_Coarsen: %d modifications in pass %d.\n",theMod->idx,iter_mods,pass);
#endif
#ifdef DEBUG_QUALITY
    tests(false);
#endif
  }
#ifdef ADAPT_VERBOSE
  CkPrintf("[%d]ParFUM_Coarsen: %d total modifications over %d passes.\n",theMod->idx,mods,pass);
#endif
  delete[] coarsenElements;
  return mods;
}

/** Performs a sequence of refinements or coarsening as is needed
 * to achieve the target areas for elements
 */
void FEM_Adapt_Algs::FEM_AdaptMesh(int qm, int method, double factor, 
				   double *sizes)
{
  MPI_Comm comm=(MPI_Comm)FEM_chunk::get("FEM_Update_mesh")->defaultComm;
  MPI_Barrier(comm);
#ifdef ADAPT_VERBOSE
  CkPrintf("[%d]BEGIN: FEM_AdaptMesh...\n",theMod->idx);
#endif
#ifdef DEBUG_QUALITY
  tests(true);
#endif
  SetMeshSize(method, factor, sizes);
  MPI_Barrier(comm);
  GradateMesh(GRADATION);
  MPI_Barrier(comm);
  (void)Refine(qm, method, factor, sizes);
  MPI_Barrier(comm);
  GradateMesh(GRADATION);
  MPI_Barrier(comm);
  (void)Coarsen(qm, method, factor, sizes);
#ifdef DEBUG_QUALITY
  MPI_Barrier(comm);
#endif
  FEM_Repair(qm);
  MPI_Barrier(comm);
#ifdef DEBUG_QUALITY
  tests(true);
  MPI_Barrier(comm);
#endif
#ifdef ADAPT_VERBOSE
  CkPrintf("[%d]...END: FEM_AdaptMesh.\n",theMod->idx);
#endif
}

/** Smooth the mesh using method according to some quality measure qm
 */
void FEM_Adapt_Algs::FEM_Smooth(int qm, int method)
{
  CkPrintf("WARNING: ParFUM_Smooth: Not yet implemented.\n");
}

/** FEM_Mesh_smooth
 *	Inputs	: meshP - a pointer to the FEM_Mesh object to smooth
 *		: nodes - an array of local node numbers to be smoothed.  Send
 *			  NULL pointer to smooth all nodes.
 *		: nNodes - the size of the nodes array
 *		: attrNo - the attribute number where the coords are registered
 *	Shifts nodes around to improve mesh quality.  FEM_BOUNDARY attribute
 *	and interpolator function must be registered by user to maintain 
 *	boundary information.
 */
void  FEM_Adapt_Algs::FEM_mesh_smooth(FEM_Mesh *meshP, int *nodes, int nNodes, int attrNo) {
  vector2d newPos, *coords, *ghostCoords;
  int idx, nNod, nGn, gIdxN, *boundVals, nodesInChunk, mesh;
  int *adjnodes;

  mesh=FEM_Mesh_default_read();
  nodesInChunk = FEM_Mesh_get_length(mesh,FEM_NODE);
  nGn = FEM_Mesh_get_length(mesh, FEM_GHOST + FEM_NODE);
  
  boundVals = new int[nodesInChunk];
  coords = new vector2d[nodesInChunk+nGn];

  FEM_Mesh_data(mesh, FEM_NODE, FEM_BOUNDARY, (int*) boundVals, 0, nodesInChunk, FEM_INT, 1);    

  FEM_Mesh_data(mesh, FEM_NODE, attrNo, (double*)coords, 0, nodesInChunk, FEM_DOUBLE, 2);

  IDXL_Layout_t coord_layout = IDXL_Layout_create(IDXL_DOUBLE, 2);
  FEM_Update_ghost_field(coord_layout,-1, coords); 
  ghostCoords = &(coords[nodesInChunk]);
  for (int i=0; i<nNodes; i++)
  {
    if (nodes==NULL) idx=i;
    else idx=nodes[i];
    newPos.x=0;
    newPos.y=0;
    CkAssert(idx<nodesInChunk);  
    if (FEM_is_valid(mesh, FEM_NODE, idx) && boundVals[idx]==0) //node must be internal
    {
      meshP->n2n_getAll(idx, adjnodes, nNod);
      for (int j=0; j<nNod; j++) { //for all adjacent nodes, find coords
	if (adjnodes[j]<-1) {
	  gIdxN = FEM_From_ghost_index(adjnodes[j]);
	  newPos.x += ghostCoords[gIdxN].x;
	  newPos.y += ghostCoords[gIdxN].y;
	}
	else {
	  newPos.x += coords[adjnodes[j]].x;
	  newPos.y += coords[adjnodes[j]].y;
	}     
      }
      newPos.x/=nNod;
      newPos.y/=nNod;
      FEM_set_entity_coord2(mesh, FEM_NODE, idx, newPos.x, newPos.y);
      delete [] adjnodes;
    }
  }

  delete [] coords;
  delete [] boundVals;
}

/** Elements with a bad quality metric are either flipped or coarsened to 
 * newer elements which are of better quality
 */
void FEM_Adapt_Algs::FEM_Repair(int qm)
{
  double avgQual = 0.0, minQual = getAreaQuality(0);
  int numBadElems = 0;
  int elemConn[3];
#ifdef ADAPT_VERBOSE
  CkPrintf("[%d]WARNING: ParFUM_Repair: Under construction.\n",theMod->idx);
  numElements = theMesh->elem[0].size();
  for (int i=0; i<numElements; i++) { 
    if (theMesh->elem[0].is_valid(i)) {
      double qFactor=getAreaQuality(i);
      avgQual += qFactor;
      if (qFactor <  QUALITY_MIN) {
	numBadElems++;
	if (qFactor < minQual) minQual = qFactor;
      }
    }
  }
  avgQual /= numElements;
  CkPrintf("BEFORE FEM_Repair: Average Element Quality = %2.6f, Min = %2.6f (1.0 is perfect)\n", avgQual, minQual);
  //CkPrintf("BEFORE FEM_Repair: Average Element Quality = %2.6f, Min = %2.6f (1.0 is perfect)\n  %d out of %d elements were below the minimum quality tolerance of %2.6f\n", avgQual, minQual, numBadElems, numElements, QUALITY_MIN);
#endif

  int elemWidth = theMesh->elem[0].getConn().width();
  int changes=1, totalChanges=0;
  int count=0;
  while (changes!=0 && count<4) {
    count++;
    changes = 0;
    numElements = theMesh->elem[0].size();
    for (int i=0; i<numElements; i++) { 
      if (theMesh->elem[0].is_valid(i)) {
	double qFactor=getAreaQuality(i);
	if (qFactor <  0.75*QUALITY_MIN) {
	  int elId = i;
	  theMesh->e2n_getAll(elId, elemConn);
	  bool notclear = false;
	  for(int j=0; j<elemWidth; j++) {
	    int nd = elemConn[j];
	    if(nd>=0) notclear = theMod->fmLockN[nd].haslocks();
	    //else cannot look up because connectivity might be screwed up,
	    //I am hoping that at least one of the nodes will be local
	  }
	  if(notclear) continue;
	  int n1=elemConn[0], n2=elemConn[1];
	  //too bad.. should not decide without locking!!
	  //these values might change by the time lock is acquired
	  double len1 = length(elemConn[0],elemConn[1]);
	  double len2 = length(elemConn[1],elemConn[2]);
	  double len3 = length(elemConn[2],elemConn[0]);
	  if(len1<-1.0 || len2<-1.0 || len3<-1.0) continue;
	  double avglen=(len1+len2+len3)/3.0;
	  int maxn1=0, maxn2=1;
	  int minn1=0, minn2=1;
	  double maxlen=len1;
	  int maxed=0;
	  if(len2>maxlen) {
	    maxlen = len2;
	    maxn1 = 1;
	    maxn2 = 2;
	    maxed = 1;
	  }
	  if(len3>maxlen) {
	    maxlen = len3;
	    maxn1 = 2;
	    maxn2 = 0;
	    maxed = 2;
	  }
	  double minlen = len1;
	  if(len2<minlen) {
	    minlen = len2;
	    minn1 = 1;
	    minn2 = 2;
	  }
	  if(len3<minlen) {
	    minlen = len3;
	    minn1 = 2;
	    minn2 = 0;
	  }
	  double otherlen = 3*avglen - maxlen - minlen;
	  if (maxlen > 0.95*(minlen+otherlen)) { // refine
	    //decide if this should be a bisect or flip, 
	    //depends on if the longest edge on this element is also the longest on the
	    //element sharing this edge, if not, flip it
	    int nbrEl = theMesh->e2e_getNbr(i,maxed);
	    double len4=0.0, len5=0.0;
	    if(nbrEl!=-1) {
	      int con1[3];
	      theMesh->e2n_getAll(nbrEl,con1);
	      int nbrnode=-1;
	      for(int j=0; j<3; j++) {
		if(con1[j]!=elemConn[maxn1] && con1[j]!=elemConn[maxn2]) {
		  nbrnode = con1[j];
		  break;
		}
	      }
	      len4 = length(elemConn[maxn1],nbrnode);
	      len5 = length(elemConn[maxn2],nbrnode);
	    }
	    int success = -1;
	    if(len4>maxlen || len5>maxlen) {
#ifdef DEBUG_FLIP
	      success = theAdaptor->edge_flip(elemConn[maxn1], elemConn[maxn2]);
#endif
	    }
	    else {
	      success = theAdaptor->edge_bisect(elemConn[maxn1], elemConn[maxn2]);
	    }
	    if (success >= 0) {
	      //CkPrintf("Refined bad element!\n");
	      changes++;
	    }
	  }
	  else if (minlen < 0.10*(maxlen+otherlen)) { // coarsen
	    int success = theAdaptor->edge_contraction(elemConn[minn1], elemConn[minn2]);
	    if (success >= 0) { 
	      //CkPrintf("Coarsened bad element!\n");
	      changes++;
	    }
	  }
	  else {
	    //CkPrintf("Leaving one bad element alone...\n");
	  }
	}
      }
    }
    totalChanges += changes;
  }

#ifdef ADAPT_VERBOSE
  numElements = theMesh->elem[0].size();
  numBadElems = 0;
  avgQual = 0.0;
  minQual = getAreaQuality(0);
  for (int i=0; i<numElements; i++) { 
    if (theMesh->elem[0].is_valid(i)) {
      double qFactor=getAreaQuality(i);
      avgQual += qFactor;
      if (qFactor <  QUALITY_MIN) {
	numBadElems++;
	if (qFactor < minQual) minQual = qFactor;
      }
    }
  }
  avgQual /= numElements;
  CkPrintf("[%d]AFTER FEM_Repair: Average Element Quality = %2.6f, Min = %2.6f (1.0 is perfect) No. of repairs %d\n",theMod->idx,avgQual, minQual,totalChanges);
#ifdef DEBUG_QUALITY
    tests(false);
#endif
  //  CkPrintf("AFTER FEM_Repair: Average Element Quality = %2.6f, Min = %2.6f (1.0 is perfect)\n  %d out of %d elements were below the minimum quality tolerance of %2.6f\n", avgQual, minQual, numBadElems, numElements, QUALITY_MIN);
#endif
}

/** Remesh entire mesh according to quality measure qm
 * if method = 0, set entire mesh size to factor
 * if method = 1, keep regional mesh sizes, and scale by factor
 * if method = 2, uses sizes to size mesh by regions 
 */
void FEM_Adapt_Algs::FEM_Remesh(int qm, int method, double factor, 
				double *sizes)
{
  CkPrintf("WARNING: ParFUM_Remesh: Under construction.\n");
}



/** For each element, set its size to its average edge length
 */
void FEM_Adapt_Algs::SetReferenceMesh()
{
  // TODO: do we need to run this loop for element types other than 0?
  double avgLength = 0.0;
  int width = theMesh->elem[0].getConn().width();
  int numElements = theMesh->elem[0].size();
  int elemConn[3];
  
  for (int i=0; i<numElements; ++i, avgLength=0) {
    if(theMesh->elem[0].is_valid(i)) {
      theMesh->e2n_getAll(i, elemConn);
      for (int j=0; j<width-1; ++j) {
	avgLength += length(elemConn[j], elemConn[j+1]);
      }
      avgLength += length(elemConn[0], elemConn[width-1]);
      avgLength /= width;
      theMesh->elem[0].setMeshSizing(i, avgLength);      
    }
  }
}

/** Resize mesh elements to avoid jumps in element size
    Algorithm based on h-shock correction, described in
    Mesh Gradation Control, Borouchaki et al
    IJNME43 1998 www.ann.jussieu.fr/~frey/publications/ijnme4398.pdf 
*/
void FEM_Adapt_Algs::GradateMesh(double smoothness)
{
    const double beta = smoothness;

    double maxShock, minShock;
    int iteration = 0, updates = 0;

    int* adjNodes, *boundNodes;
    int nadjNodes, nnodes;
    int meshNum = FEM_Mesh_default_read();
    //if (smoothness < 1.0) {
      //printf("");
    //}

    nnodes = theMesh->node.size();
    boundNodes = new int[nnodes];
    FEM_Mesh_data(meshNum, FEM_NODE, FEM_BOUNDARY, 
            boundNodes, 0, nnodes, FEM_INT, 1);


    //printf("Running h-shock mesh gradation with beta=%.3f\n", beta);
    fflush(NULL);

#ifndef GRADATION_ITER_LIMIT
#define GRADATION_ITER_LIMIT    10
#endif
    
    do {
        maxShock = 0;
        minShock = 1e10;
        
        for (int node=0; node<nnodes; ++node) {
            if (boundNodes[node]!= 0 || !FEM_is_valid(meshNum, FEM_NODE, node))
               continue;
            //if (!FEM_is_valid(meshNum, FEM_NODE, node))
            //    continue;
            
            theMesh->n2n_getAll(node, adjNodes, nadjNodes);
            for (int adjNode=0; adjNode<nadjNodes; ++adjNode) {
                double edgelen = length(node, adjNodes[adjNode]);
                
                // get adjacent elemnents and their sizes
                int e1, e2;
                theMesh->get2ElementsOnEdge(node, adjNodes[adjNode], &e1, &e2);
                if (e1 <= -1 || e2 <= -1) continue; //this will not smooth across boundaries
                
                double s1, s2;
                s1 = theMesh->elem[0].getMeshSizing(e1);
                s2 = theMesh->elem[0].getMeshSizing(e2);
                if (s1 <= 0 || s2 <= 0) continue;
                
                // h-shock=max(size ratio)^(1/edge length)
                CkAssert(s1 >= 0 && s2 >= 0 && "Bad size");
                //CkAssert(edgelen > 1e-6 && "Length 0 edge");
                CkAssert(edgelen == edgelen && "Length inf edge");

                double ratio = (s1 > s2) ? s1/s2 : s2/s1;
                CkAssert (ratio >= 1.0 && ratio == ratio && "Bad ratio");
                
                // WARNING WARNING WARNING
                // TEST ONLY, THIS IS NOT CORRECT
                if (ratio > beta) {
                    if (s1 > s2) {
                        theMesh->elem[0].setMeshSizing(e1, s1 - (s1-s2)/3);
                        theMesh->elem[0].setMeshSizing(e2, s2 + (s1-s2)/3);
                    } else {
                        theMesh->elem[0].setMeshSizing(e2, s2 - (s2-s1)/3);
                        theMesh->elem[0].setMeshSizing(e1, s1 + (s2-s1)/3);
                    }
                    updates++;
                }
                if (ratio > maxShock) maxShock = ratio;
                if (ratio < minShock) minShock = ratio;
                
                
                ////double hs = ratio;
                //double hs = pow(ratio, 1.0/edgelen);
                //
                //if (hs > maxShock) maxShock = hs;
                //if (hs < minShock) minShock = hs;

                //// if hs > beta, resize the larger elt:
                //// new size = old size / eta^2
                //// eta = (beta / h-shock)^(edge length)
                ////     = beta^(edge length) / size ratio
                //
                //if (hs > beta) {
                //    double etasq = pow(beta, edgelen) / ratio;
                //    etasq *= etasq;

                //    //if (hs > 100) {
                //    //    printf("hs: %8.5f\ns1: %8.5f\ns2: %8.5f\nedgelen: %8.5f\nratio: %8.5f\netasq: %8.5f", hs, s1, s2, edgelen, ratio, etasq);
                //    //    abort();
                //    //}
                //    
                //    if (s1 > s2) {
                //        theMesh->elem[0].setMeshSizing(e1, s1 / etasq);
                //    } else {
                //        theMesh->elem[0].setMeshSizing(e2, s2 / etasq);
                //    }
                //    updates++;
                //}
                
            }
            free(adjNodes);
        } 
        
	//printf("Finished iteration %d\n", iteration);
        //printf("Max shock:%8.3f\n", maxShock);
        //printf("Min shock:%8.3f\n", minShock);
        //printf("Target:%8.3f\n", beta);
        
    } while (maxShock > beta && ++iteration < GRADATION_ITER_LIMIT);

    //printf("%d total updates in %d iterations in GradateMesh\n", updates, iteration);
    fflush(NULL);

    delete[] boundNodes;

    return;
}

/** Set sizes on elements throughout the mesh; note: size is edge length */
void FEM_Adapt_Algs::SetMeshSize(int method, double factor, double *sizes)
{
  numNodes = theMesh->node.size();
  numElements = theMesh->elem[0].size();
  int elemConn[3];

  if (method == 0) { // set uniform sizing specified in factor

    for (int i=0; i<numElements; i++) {
      theMesh->elem[0].setMeshSizing(i, factor);
    }
    //CkPrintf("ParFUM_SetMeshSize: UNIFORM %4.6e\n", factor);
  }
  else if (method == 1) { // copy sizing from array
    for (int i=0; i<numElements; i++) {
      if (sizes[i] > 0.0) {
	theMesh->elem[0].setMeshSizing(i, sizes[i]);
      }
    }
    //CkPrintf("ParFUM_SetMeshSize: SIZES input\n");
  }
  else if (method == 2) { // calculate current sizing and scale by factor
    double avgEdgeLength = 0.0;
    int width = theMesh->elem[0].getConn().width();
    int numEdges=3;
    if (dim==3) numEdges=6;
    for (int i=0; i<numElements; i++) {
      if(theMesh->elem[0].is_valid(i)) {
	theMesh->e2n_getAll(i, elemConn);
	for (int j=0; j<width-1; j++) {
	  for (int k=j+1; k<width; k++) {
	    avgEdgeLength += length(elemConn[j], elemConn[k]);
	  }
	}
	avgEdgeLength += length(elemConn[0], elemConn[width-1]);
	avgEdgeLength /= (double)numEdges;
	theMesh->elem[0].setMeshSizing(i, factor*avgEdgeLength);
      }
    }
    //CkPrintf("ParFUM_SetMeshSize: CALCULATED & SCALED \n");
  }
  else if (method == 3) { // scale existing sizes by array sizes
    for (int i=0; i<numElements; i++) {
      if (sizes[i] > 0.0) {
	theMesh->elem[0].setMeshSizing(i, sizes[i]*theMesh->elem[0].getMeshSizing(i));
      }
    }
  }
  else if (method == 4) { // scale existing sizes by factor
    for (int i=0; i<numElements; i++) {
      theMesh->elem[0].setMeshSizing(i, factor*theMesh->elem[0].getMeshSizing(i));
    }
  }
  else if (method == 5) { // mesh sizing has been set independently; use as is
    //CkPrintf("ParFUM_SetMeshSize: USE EXISTING SIZES \n");
  }
  //  CkPrintf("Current mesh sizing: ");
  //for (int i=0; i<numElements; i++) {
  //CkPrintf("%4.6e ", theMesh->elem[0].getMeshSizing(i));
  //}
}

/** Insert this element to the refine or coarsen heap
 */
void FEM_Adapt_Algs::Insert(int eIdx, double len, int cFlag)
{
  int i;
  if (cFlag) {
    i = ++coarsenHeapSize; 
    while ((coarsenElements[i/2].len>=len) && (i != 1)) {
      coarsenElements[i].len=coarsenElements[i/2].len;
      coarsenElements[i].elID=coarsenElements[i/2].elID;
      i/=2;
    }
    coarsenElements[i].elID=eIdx;
    coarsenElements[i].len=len; 
  }
  else {
    i = ++refineHeapSize; 
    while ((refineElements[i/2].len>=len) && (i != 1)) {
      refineElements[i].len=refineElements[i/2].len;
      refineElements[i].elID=refineElements[i/2].elID;
      i/=2;
    }
    refineElements[i].elID=eIdx;
    refineElements[i].len=len; 
  }
}

/** Removes and returns the minimum element from the refine/coarsen heap
 */
int FEM_Adapt_Algs::Delete_Min(int cflag)
{
  int Child, i, Min_ID; 
  if (cflag) {
    Min_ID=coarsenElements[1].elID;
    for (i=1; i*2 <= coarsenHeapSize-1; i=Child) { // Find smaller child
      Child = i*2; // child is left child  
      if (Child != coarsenHeapSize)  // right child exists
	if (coarsenElements[Child+1].len < coarsenElements[Child].len)
	  Child++; 
      // Percolate one level
      if (coarsenElements[coarsenHeapSize].len >= coarsenElements[Child].len) {
	coarsenElements[i].elID = coarsenElements[Child].elID;
	coarsenElements[i].len = coarsenElements[Child].len;
      }
      else break; 
    }
    coarsenElements[i].elID = coarsenElements[coarsenHeapSize].elID;
    coarsenElements[i].len = coarsenElements[coarsenHeapSize].len; 
    coarsenHeapSize--;
    return Min_ID; 
  }
  else {
    Min_ID=refineElements[1].elID;
    for (i=1; i*2 <= refineHeapSize-1; i=Child) { // Find smaller child
      Child = i*2;       // child is left child  
      if (Child !=refineHeapSize)  // right child exists
	if (refineElements[Child+1].len < refineElements[Child].len)
	  Child++; 
      // Percolate one level
      if (refineElements[refineHeapSize].len >= refineElements[Child].len){  
	refineElements[i].elID = refineElements[Child].elID;   
	refineElements[i].len = refineElements[Child].len;
      }
      else break; 
    }
    refineElements[i].elID = refineElements[refineHeapSize].elID;
    refineElements[i].len = refineElements[refineHeapSize].len; 
    refineHeapSize--;
    return Min_ID; 
  }
}



// =====================  BEGIN refine_element_leb ========================= 
/** Initiate instance of Longest Edge Bisection on element e.  Propagates
    throughout the mesh to maintain the requirement that only longest edges
    are bisected; returns the new node index that it created
    Given an element e, if e's longest edge f is also the longest edge
    of e's neighbor across f, g, split f by adding a new node in the 
    center of f, and splitting both e and g into two elements.  If g
    does not have f as it's longest edge, recursively call refine_element_leb 
    on g, and start over. 
*/ 
int FEM_Adapt_Algs::refine_element_leb(int e) {
  int fixNode, otherNode, opNode, longEdge, nbr; 
  double eLens[3], longEdgeLen = 0.0;
  int elemConn[3];

  if(e==-1) {
    return -1;
  }

  theMesh->e2n_getAll(e, elemConn);
  eLens[0] = length(elemConn[0], elemConn[1]);
  eLens[1] = length(elemConn[1], elemConn[2]);
  eLens[2] = length(elemConn[2], elemConn[0]);
  for (int i=0; i<3; i++) {
    if (eLens[i] > longEdgeLen) {
      longEdgeLen = eLens[i];
      longEdge = i;
      fixNode = elemConn[i];
      otherNode = elemConn[(i+1)%3];
      opNode = elemConn[(i+2)%3];
    }
  }

  nbr = theMesh->e2e_getNbr(e, longEdge);
  if (nbr == -1) // e's longEdge is on physical boundary
    return theAdaptor->edge_bisect(fixNode, otherNode);
  int nbrOpNode = theAdaptor->e2n_getNot(nbr, fixNode, otherNode);
  double fixEdgeLen = length(fixNode, nbrOpNode);
  double otherEdgeLen = length(otherNode, nbrOpNode);
  if ((fixEdgeLen > longEdgeLen) || (otherEdgeLen > longEdgeLen)) { 
    // longEdge is not nbr's longest edge
    int newNode = theAdaptor->edge_bisect(fixNode, otherNode);
    if(newNode==-1) return -1;
    int propElem, propNode; // get the element to propagate on
    if (fixEdgeLen > otherEdgeLen) {
      propElem = theAdaptor->findElementWithNodes(newNode, fixNode, nbrOpNode);
      propNode = fixNode;
    }
    else {
      propElem = theAdaptor->findElementWithNodes(newNode,otherNode,nbrOpNode);
      propNode = otherNode;
    }

    //if propElem is ghost, then it's propagating to neighbor, otherwise not
    if(!FEM_Is_ghost_index(propElem)) {
      refine_flip_element_leb(propElem,propNode,newNode,nbrOpNode,longEdgeLen);
    }
    else {
      int localChk, nbrChk;
      localChk = theMod->getfmUtil()->getIdx();
      nbrChk = theMod->getfmUtil()->getRemoteIdx(theMesh,propElem,0);
      int propNodeT = theAdaptor->getSharedNodeIdxl(propNode, nbrChk);
      int newNodeT = theAdaptor->getSharedNodeIdxl(newNode, nbrChk);
      int nbrghost = (nbrOpNode>=0)?0:1;
      int nbrOpNodeT = (nbrOpNode>=0)?(theAdaptor->getSharedNodeIdxl(nbrOpNode, nbrChk)):(theAdaptor->getGhostNodeIdxl(nbrOpNode, nbrChk));
      int propElemT = theAdaptor->getGhostElementIdxl(propElem, nbrChk);
      meshMod[nbrChk].refine_flip_element_leb(localChk, propElemT, propNodeT, newNodeT,nbrOpNodeT,nbrghost,longEdgeLen);
    }
    return newNode;
  }
  else return theAdaptor->edge_bisect(fixNode, otherNode); // longEdge on nbr
}

/** Continue longest edge bisection on element 'e', then flip edge (n1,n2);
    If the length of the new edge (p, newnode) is greater than 'le' (longest edge),
    recursively call this function on (newElem, p, n1, newNode, le)
*/
void FEM_Adapt_Algs::refine_flip_element_leb(int e, int p, int n1, int n2, 
					     double le) 
{
  int newNode = refine_element_leb(e);
  if(newNode == -1) return;
  (void) theAdaptor->edge_flip(n1, n2);
  if (length(p, newNode) > le) {
    int localChk = theMod->getfmUtil()->getIdx();
    int newElem = theAdaptor->findElementWithNodes(newNode, n1, p);
    refine_flip_element_leb(newElem, p, n1, newNode, le);
  }
}
// ========================  END refine_element_leb ========================



/** The region is defined by (xmax, xmin, ymax, ymin) and the target area is given by targetA
 */
int FEM_Adapt_Algs::simple_refine(double targetA, double xmin, double ymin, double xmax, double ymax) {
  bool adapted = true;
  refineElements = refineStack = NULL;
  refineTop = refineHeapSize = 0;
  int elemConn[3];
  double coordsn1[2], coordsn2[2], coordsn3[2];
  int elemWidth=3;

  while(adapted) {
    adapted = false;
    int noEle = theMesh->elem[0].size();
    refineElements = new elemHeap[noEle+1];
    for(int i=0; i<noEle; i++) {
      if(theMesh->elem[0].is_valid(i)) {
	theMesh->e2n_getAll(i,elemConn,0);
	bool notclear = false;
	for(int j=0; j<elemWidth; j++) {
	  int nd = elemConn[j];
	  if(nd>=0) notclear = theMod->fmLockN[nd].haslocks();
	  //else cannot look up because connectivity might be screwed up,
	  //I am hoping that at least one of the nodes will be local
	}
	if(notclear) continue;
	getCoord(elemConn[0], coordsn1);
	getCoord(elemConn[1], coordsn2);
	getCoord(elemConn[2], coordsn3);
	if(coordsn1[0]<-1.0 ||coordsn1[1]<-1.0 || coordsn2[0]<-1.0 ||coordsn2[1]<-1.0 || coordsn3[0]<-1.0 ||coordsn3[1]<-1.0) continue;
	//do a refinement only if it has any node within the refinement box
	if((coordsn1[0]<xmax && coordsn1[0]>xmin && coordsn1[1]<ymax && coordsn1[1]>ymin) || (coordsn2[0]<xmax && coordsn2[0]>xmin && coordsn2[1]<ymax && coordsn2[1]>ymin) || (coordsn3[0]<xmax && coordsn3[0]>xmin && coordsn3[1]<ymax && coordsn3[1]>ymin)) {
	  Insert(i,-getArea(coordsn1, coordsn2, coordsn3),0); 
	  //negative area to trick the minheap to become a maxheap!
	}
	else {
	  Insert(i,-MINAREA,0); //make it believe that this triangle does not need refinement
	}
      } else {
	Insert(i,-MINAREA,0);
      }
    }

    //bad bubblesort
    /*for(int i=0; i<noEle; i++) {
      for(int j=i+1; j<noEle; j++) {
	if(areas[j] > areas[i]) {
	  double tmp = areas[j];
	  areas[j] = areas[i];
	  areas[i] = tmp;
	  int t = map1[j];
	  map1[j] = map1[i];
	  map1[i] = t;
	}
      }
      }*/
    
    /*for(int i=0; i<noEle; i++) {
      if(theMesh->elem[0].is_valid(map1[i])) {
	if(areas[i] > targetA) {
	  int ret = refine_element_leb(map1[i]);
	  if(ret!=-1) {
	    adapted = true;
	  }
	}
	}
      //if(adapted) break;
      }*/
    while(refineHeapSize>0) {
      int tmp = Delete_Min(0);
      if(tmp!=-1 && theMesh->elem[0].is_valid(tmp)) {
	//since we are reusing element indices, so do not trust earlier areas
	theMesh->e2n_getAll(tmp,elemConn,0);
	bool notclear = false;
	for(int j=0; j<3; j++) {
	  int nd = elemConn[j];
	  if(nd>=0) notclear = theMod->fmLockN[nd].haslocks();
	  //else cannot look up because connectivity might be screwed up,
	  //I am hoping that at least one of the nodes will be local
	}
	if(notclear) continue;
	getCoord(elemConn[0], coordsn1);
	getCoord(elemConn[1], coordsn2);
	getCoord(elemConn[2], coordsn3);
	if(coordsn1[0]<-1.0 ||coordsn1[1]<-1.0 || coordsn2[0]<-1.0 ||coordsn2[1]<-1.0 || coordsn3[0]<-1.0 ||coordsn3[1]<-1.0) continue;
	if(getArea(coordsn1, coordsn2, coordsn3) > targetA) {
	  double len1 = length(coordsn1,coordsn2);
	  double len2 = length(coordsn2,coordsn3);
	  double len3 = length(coordsn3,coordsn1);
	  int n1=elemConn[0], n2=elemConn[1];
	  double max=len1;
	  if(len2>max) {
	    max = len2; n1 = elemConn[1]; n2 = elemConn[2];
	  }
	  if(len3>max) {
	    max = len3; n1 = elemConn[2]; n2 = elemConn[0];
	  }
	  //int ret = refine_element_leb(tmp);
	  int ret = theAdaptor->edge_bisect(n1,n2);
	  if(ret!=-1) adapted = true;
	}
      }
      //if(adapted) break;
    }
    delete[] refineElements;
    //if(adapted) break;
  }

  if(adapted) return -1;
  else return 1;
}

/** The region is defined by (xmax, xmin, ymax, ymin) and the target area is given by targetA
 */
int FEM_Adapt_Algs::simple_coarsen(double targetA, double xmin, double ymin, double xmax, double ymax) {
  int noEle = theMesh->elem[0].size();
  int *shortestEdge = (int *)malloc(2*sizeof(int));
  bool adapted = true;
  coarsenElements = NULL;
  coarsenHeapSize = 0;
  int elemConn[3];
  double coordsn1[2], coordsn2[2], coordsn3[2];
  int elemWidth=3;

  while(adapted) {
    adapted = false;
    coarsenElements = new elemHeap[noEle+1];
    coarsenElements[0].len=-2.0;
    coarsenElements[0].elID=-1;
    for(int i=0; i<noEle; i++) {
      if(theMesh->elem[0].is_valid(i)) {
	theMesh->e2n_getAll(i,elemConn,0);
	bool notclear = false;
	for(int j=0; j<elemWidth; j++) {
	  int nd = elemConn[j];
	  if(nd>=0) notclear = theMod->fmLockN[nd].haslocks();
	  //else cannot look up because connectivity might be screwed up,
	  //I am hoping that at least one of the nodes will be local
	}
	if(notclear) continue;
	getCoord(elemConn[0], coordsn1);
	getCoord(elemConn[1], coordsn2);
	getCoord(elemConn[2], coordsn3);
	if(coordsn1[0]<-1.0 ||coordsn1[1]<-1.0 || coordsn2[0]<-1.0 ||coordsn2[1]<-1.0 || coordsn3[0]<-1.0 ||coordsn3[1]<-1.0) continue;
	//do a coarsening only if it has any node within the coarsen box
	if((coordsn1[0]<xmax && coordsn1[0]>xmin && coordsn1[1]<ymax && coordsn1[1]>ymin) || (coordsn2[0]<xmax && coordsn2[0]>xmin && coordsn2[1]<ymax && coordsn2[1]>ymin) || (coordsn3[0]<xmax && coordsn3[0]>xmin && coordsn3[1]<ymax && coordsn3[1]>ymin)) {
	  Insert(i,getArea(coordsn1, coordsn2, coordsn3),1);
	} 
	else {
	  Insert(i,MAXAREA,1); //make it believe that this triangle is big enough
	}
      } else {
	Insert(i,MAXAREA,1);
      }
    }
    
    /*for(int i=0; i<noEle; i++) {
      for(int j=i+1; j<noEle; j++) {
	if(areas[j] < areas[i]) {
	  double tmp = areas[j];
	  areas[j] = areas[i];
	  areas[i] = tmp;
	  int t = map1[j];
	  map1[j] = map1[i];
	  map1[i] = t;
	}
      }
      }*/
    
    /*for(int i=0; i<noEle; i++) {
      if(theMesh->elem[0].is_valid(map1[i])) {
	if(areas[i] < targetA) {
	  //find the nodes along the smallest edge & coarsen the edge
	  theMesh->e2n_getAll(map1[i],elemConn,0);
	  getShortestEdge(elemConn[0], elemConn[1], elemConn[2], shortestEdge);
	  int ret = theAdaptor->edge_contraction(shortestEdge[0], shortestEdge[1]);
	  if(ret != -1) adapted = true;
	}
      }
      //if(adapted) break;
      }*/
    while(coarsenHeapSize>0) {
      int tmp = Delete_Min(1);
      if(tmp!=-1 && theMesh->elem[0].is_valid(tmp)) {
	theMesh->e2n_getAll(tmp,elemConn,0);
	bool notclear = false;
	for(int j=0; j<3; j++) {
	  int nd = elemConn[j];
	  if(nd>=0) notclear = theMod->fmLockN[nd].haslocks();
	  //else cannot look up because connectivity might be screwed up,
	  //I am hoping that at least one of the nodes will be local
	}
	if(notclear) continue;
	getCoord(elemConn[0], coordsn1);
	getCoord(elemConn[1], coordsn2);
	getCoord(elemConn[2], coordsn3);
	if(coordsn1[0]<-1.0 ||coordsn1[1]<-1.0 || coordsn2[0]<-1.0 ||coordsn2[1]<-1.0 || coordsn3[0]<-1.0 ||coordsn3[1]<-1.0) continue;
	if(getArea(coordsn1, coordsn2, coordsn3) < targetA) {
	  getShortestEdge(elemConn[0], elemConn[1], elemConn[2], shortestEdge);
	  int ret = theAdaptor->edge_contraction(shortestEdge[0], shortestEdge[1]);
	  if(ret != -1) adapted = true;
	}
      }
      //if(adapted) break;
    }
    delete[] coarsenElements;
    //if(adapted) break;
  }
  free(shortestEdge);

  if(adapted) return -1;
  else return 1;
}



//HELPER functions
double FEM_Adapt_Algs::length(int n1, int n2)
{
  double coordsn1[2], coordsn2[2];
  getCoord(n1, coordsn1);
  getCoord(n2, coordsn2);

  double ret = length(coordsn1, coordsn2);
  return ret;
}

double FEM_Adapt_Algs::length(double *n1_coord, double *n2_coord) { 
  double d, ds_sum=0.0;

  for (int i=0; i<dim; i++) {
    if(n1_coord[i]<-1.0 || n2_coord[i]<-1.0) return -2.0;
    d = n1_coord[i] - n2_coord[i];
    ds_sum += d*d;
  }
  return (sqrt(ds_sum));
}

double FEM_Adapt_Algs::getArea(int n1, int n2, int n3)
{
  double coordsn1[2], coordsn2[2], coordsn3[2];
  getCoord(n1, coordsn1);
  getCoord(n2, coordsn2);
  getCoord(n3, coordsn3);

  double ret = getArea(coordsn1, coordsn2, coordsn3);
  return ret;
}

double FEM_Adapt_Algs::getArea(double *n1_coord, double *n2_coord, double *n3_coord) {
  double area=0.0;
  double aLen, bLen, cLen, sLen, d, ds_sum;

  ds_sum = 0.0;
  for (int i=0; i<2; i++) {
    d = n1_coord[i] - n2_coord[i];
    ds_sum += d*d;
  }
  aLen = sqrt(ds_sum);
  ds_sum = 0.0;
  for (int i=0; i<2; i++) {
    d = n2_coord[i] - n3_coord[i];
    ds_sum += d*d;
  }
  bLen = sqrt(ds_sum);
  ds_sum = 0.0;
  for (int i=0; i<2; i++) {
    d = n3_coord[i] - n1_coord[i];
    ds_sum += d*d;
  }
  cLen = sqrt(ds_sum);
  sLen = (aLen+bLen+cLen)/2;
  if(sLen-aLen < 0) return 0.0;
  else if(sLen-bLen < 0) return 0.0;
  else if(sLen-cLen < 0) return 0.0; //area too small to note
  return (sqrt(sLen*(sLen-aLen)*(sLen-bLen)*(sLen-cLen)));
}

double FEM_Adapt_Algs::getSignedArea(int n1, int n2, int n3) {
  double coordsn1[2], coordsn2[2], coordsn3[2];
  getCoord(n1, coordsn1);
  getCoord(n2, coordsn2);
  getCoord(n3, coordsn3);

  double ret = getSignedArea(coordsn1, coordsn2, coordsn3);
  return ret;
}

double FEM_Adapt_Algs::getSignedArea(double *n1_coord, double *n2_coord, double *n3_coord) {
  double area=0.0;
  double vec1_x, vec1_y, vec2_x, vec2_y;

  vec1_x = n1_coord[0] - n2_coord[0];
  vec1_y = n1_coord[1] - n2_coord[1];
  vec2_x = n3_coord[0] - n2_coord[0];
  vec2_y = n3_coord[1] - n2_coord[1];

  area = vec1_x*vec2_y - vec2_x*vec1_y;
  return area;
}

int FEM_Adapt_Algs::getCoord(int n1, double *crds) {
  if(!FEM_Is_ghost_index(n1)) {
    FEM_Mesh_dataP(theMesh, FEM_NODE, coord_attr, (void *)crds, n1, 1, FEM_DOUBLE, dim);
  }
  else {
    int numchunks;
    IDXL_Share **chunks1;
    theMod->fmUtil->getChunkNos(0,n1,&numchunks,&chunks1);
    int index = theMod->idx;
    for(int j=0; j<numchunks; j++) {
      int chk = chunks1[j]->chk;
      if(chk==index) continue;
      int ghostidx = theMod->fmUtil->exists_in_IDXL(theMesh,n1,chk,2);
      double2Msg *d = meshMod[chk].getRemoteCoord(index,ghostidx);
      crds[0] = d->i;
      crds[1] = d->j;
      for(int j=0; j<numchunks; j++) {
	delete chunks1[j];
      }
      if(numchunks != 0) free(chunks1);
      delete d;
      break;
    }
  }
  return 1;
}

int FEM_Adapt_Algs::getShortestEdge(int n1, int n2, int n3, int* shortestEdge) {
  //note that getCoord might be a remote call, which means 
  //it might not have the same value in the mem
  //location if the memory is reused by someone else meanwhile!
  double coordsn1[2], coordsn2[2], coordsn3[2];
  getCoord(n1, coordsn1);
  getCoord(n2, coordsn2);
  getCoord(n3, coordsn3);

  double aLen = length(coordsn1, coordsn2);
  int shortest = 0;

  double bLen = length(coordsn2, coordsn3);
  if(bLen < aLen) shortest = 1;

  double cLen = length(coordsn3, coordsn1);
  if((cLen < aLen) && (cLen < bLen)) shortest = 2;

  if(shortest==0) {
    shortestEdge[0] = n1;
    shortestEdge[1] = n2;
  }
  else if (shortest==1) {
    shortestEdge[0] = n2;
    shortestEdge[1] = n3;
  }
  else {
    shortestEdge[0] = n3;
    shortestEdge[1] = n1;
  }
  return 1;
}



/** The quality metric is proportional to the area of the triangle upon
    the sum of squares of the lengths of each of the 3 sides
*/
double FEM_Adapt_Algs::getAreaQuality(int elem)
{
  double f, q, len[3];
  int n[3];
  double currentArea;
  double coordsn1[2], coordsn2[2], coordsn3[2];

  theMesh->e2n_getAll(elem, n);
  getCoord(n[0], coordsn1);
  getCoord(n[1], coordsn2);
  getCoord(n[2], coordsn3);

  currentArea = getArea(coordsn1, coordsn2, coordsn3);

  len[0] = length(coordsn1, coordsn2);
  len[1] = length(coordsn2, coordsn3);
  len[2] = length(coordsn3, coordsn1);
  f = 4.0*sqrt(3.0); //proportionality constant
  q = (f*currentArea)/(len[0]*len[0]+len[1]*len[1]+len[2]*len[2]);
  return q;
}

/** The quality metric used is the maxEdgeLength / shortestAltitude of the triangle
    If this metric is less than 100 and the area is greater than sliverArea, 
    the quality is supposed to be good
*/
void FEM_Adapt_Algs::ensureQuality(int n1, int n2, int n3) {
  double coordsn1[2], coordsn2[2], coordsn3[2];
  getCoord(n1, coordsn1);
  getCoord(n2, coordsn2);
  getCoord(n3, coordsn3);
  double area = getSignedArea(coordsn1, coordsn2, coordsn3);
  double len1 = length(coordsn1, coordsn2);
  double len2 = length(coordsn2, coordsn3);
  double len3 = length(coordsn3, coordsn1);
  double max = len1;
  if(len2>max) max = len2;
  if(len3>max) max = len3;
  //shortest edge
  double min = len1;
  if(len2<min) min = len2;
  if(len3<min) min = len3;
  double shortest_al = fabs(area/max);
  double largestR = max/shortest_al;
  CkAssert(largestR<=100.0 && -area > SLIVERAREA);
}

/** Verify the quality of the two new elements that will be created by flip
 */
bool FEM_Adapt_Algs::controlQualityF(int n1, int n2, int n3, int n4) {
  //n1 or n2 will be replaced by n4
  double coordsn1[2], coordsn2[2], coordsn3[2], coordsn4[2];
  if(n4==-1) return false;
  getCoord(n1, coordsn1);
  getCoord(n2, coordsn2);
  getCoord(n3, coordsn3);
  getCoord(n4, coordsn4);
  bool flag = false;
  if(!flag) flag = controlQualityC(coordsn3,coordsn1,coordsn2,coordsn4);
  if(!flag) flag = controlQualityC(coordsn4,coordsn2,coordsn1,coordsn3);
  return flag;
}

/** Verify the quality of the four new elements that will be created by flip
 */
bool FEM_Adapt_Algs::controlQualityR(int n1, int n2, int n3, int n4) {
  //n1 or n2 will be replaced by n4
  double coordsn1[2], coordsn2[2], coordsn3[2], coordsn4[2], coordsn5[2];
  if(n4==-1) return false;
  getCoord(n1, coordsn1);
  getCoord(n2, coordsn2);
  getCoord(n3, coordsn3);
  getCoord(n4, coordsn4);
  coordsn5[0] = (coordsn1[0]+coordsn2[0])*0.5;
  coordsn5[1] = (coordsn1[1]+coordsn2[1])*0.5;
  bool flag = false;
  if(!flag) flag = theMod->fmAdaptAlgs->controlQualityR(coordsn1,coordsn3,coordsn5);
  if(!flag) flag = theMod->fmAdaptAlgs->controlQualityR(coordsn2,coordsn3,coordsn5);
  if(!flag) flag = theMod->fmAdaptAlgs->controlQualityR(coordsn1,coordsn4,coordsn5);
  if(!flag) flag = theMod->fmAdaptAlgs->controlQualityR(coordsn2,coordsn4,coordsn5);
  return flag;
}

bool FEM_Adapt_Algs::controlQualityR(double *coordsn1, double *coordsn2, double *coordsn3) {
  double area = getArea(coordsn1, coordsn2, coordsn3);
  //do some quality preservation
  double len1 = length(coordsn1, coordsn2);
  double len2 = length(coordsn2, coordsn3);
  double len3 = length(coordsn3, coordsn1);
  //longest edge
  double max = len1;
  if(len2>max) max = len2;
  if(len3>max) max = len3;
  double shortest_al = area/max;
  double largestR = max/shortest_al;
  if(largestR>50.0) return true;
  else return false;
}

/** Computes the longestLength, shotestLength, shortestAltitude, and other
    quality metrics for the new element (n1,n2,n4)
    If the new element forms a sliver or bad quality, return true
*/
bool FEM_Adapt_Algs::controlQualityC(int n1, int n2, int n3, double *n4_coord) {
  //n3 is the node to be deleted, n4 is the new node to be added
  double coordsn1[2], coordsn2[2], coordsn3[2];
  getCoord(n1, coordsn1);
  getCoord(n2, coordsn2);
  getCoord(n3, coordsn3);
  return controlQualityC(coordsn1,coordsn2,coordsn3,n4_coord);
}

bool FEM_Adapt_Algs::controlQualityC(double *coordsn1, double *coordsn2, double *coordsn3, double *n4_coord) {
  double ret_old = getSignedArea(coordsn1, coordsn2, coordsn3);
  double ret_new = getSignedArea(coordsn1, coordsn2, n4_coord);
  //do some quality preservation
  double len1 = length(coordsn1, coordsn2);
  double len2 = length(coordsn2, n4_coord);
  double len3 = length(n4_coord, coordsn1);
  //longest edge
  double max = len1;
  if(len2>max) max = len2;
  if(len3>max) max = len3;
  //shortest edge
  double min = len1;
  if(len2<min) min = len2;
  if(len3<min) min = len3;
  double shortest_al = ret_new/max;
  double largestR = max/shortest_al;
  if(ret_old > SLIVERAREA && ret_new < -SLIVERAREA) return true; //it is a flip
  else if(ret_old < -SLIVERAREA && ret_new > SLIVERAREA) return true; //it is a flip
  else if(fabs(ret_new) < SLIVERAREA) return true; // it is a sliver
  //else if(fabs(shortest_al) < 1e-5) return true; //shortest altitude is too small
  //else if(fabs(min) < 1e-5) return true; //shortest edge is too small
  else if(fabs(largestR)>50.0) return true;
  else return false;
}



/** Used from FEM_Repair, when poor quality elements are encountered, depending
    on their maxEdgeLength, it could be better to flip it or bisect it to get
    better quality elements.
    Returns true if it should flip.
*/
bool FEM_Adapt_Algs::flipOrBisect(int elId, int n1, int n2, int maxEdgeIdx, double maxlen) {
  //return true if it should flip
  int nbrEl = theMesh->e2e_getNbr(elId,maxEdgeIdx);
  double len4=0.0, len5=0.0;
  if(nbrEl!=-1) {
    int con1[3];
    theMesh->e2n_getAll(nbrEl,con1);
    int nbrnode=-1;
    for(int j=0; j<3; j++) {
      if(con1[j]!=n1 && con1[j]!=n2) {
	nbrnode = con1[j];
	break;
      }
    }
    len4 = length(n1,nbrnode);
    len5 = length(n2,nbrnode);
  }
  if(len4>1.2*maxlen || len5>1.2*maxlen) {
    return true;
  }
  else return false;
}

void FEM_Adapt_Algs::tests(bool b=true) {
  //test the mesh for slivered triangles
  if(!b) {
    theMod->fmUtil->AreaTest(theMesh);
    return;
  }
  theMod->fmUtil->AreaTest(theMesh);
  theMod->fmUtil->StructureTest(theMesh);
  theMod->fmUtil->IdxlListTest(theMesh);
  theMod->fmUtil->residualLockTest(theMesh);
  /*for(int i=0; i<theMesh->node.size(); i++) {
    if(theMesh->node.is_valid(i)) CkPrintf("Valid -- ");
    else  CkPrintf("Invalid -- ");
    theMod->fmUtil->FEM_Print_coords(theMesh,i);
    }*/
  return;
}
